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Abstract 

We study diffusion on a substrate witli permanent traps distributed with 
critical positional correlation, modeled by their placement on the perimeters 
of a critical percolation cluster. We perform a numerical analysis of the vi- 
brational density of states and the largest eigenvalue of the equivalent scalar 
elasticity problem using the method of Arnoldi and Saad. We show that the 
critical trap correlation increases the exponent appearing in the stretched ex- 
ponential behavior of the low frequency density of states by approximately a 
factor of two as compared to the case of no correlations. A finite size scaling 
hypothesis of the largest eigenvalue is proposed and its relation to the density 
of states is given. The numerical analysis of this scaling postulate leads to 
the estimation of the stretch exponent in good agreement with the density of 
states result. 
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I. INTRODUCTION 



Understanding the behavior of a particle diffusing in the presence of traps is important 
as it captures the essence of many physical processes, including diffusion controlled reactions 
where one of the reactants is immobile, trapping of excitons, etc. [|l| Mapping of the diffusion 
problem with traps to the scalar elasticity problem also affords an insight into lattice vibra- 
tion of systems such as binary alloys with contrasting elastic constants. From a theoretical 
point of view, this problem is analogous to the ideal chain in inhomogeneous media. The 
ideal chain problem is interesting as it is the Gaussian limit of the self-avoiding walks and 
yet possesses a universality class of its own in inhomogeneous media, distinct from the usual 
random walk. 

These considerations resulted in many analytical and numerical attempts to gain under- 
standing of the problem . Typical quantities used to characterize diffusion in the presence 
of traps include Po{t), the probability that the diffusing particle returns to the starting point 
after time t. Po{t) is related to the number of distinct sites visited by the diffusing particle 
and is the Laplace transform of the vibrational density of states of the corresponding scalar 
elasticity problem, which can often be measured experimentally by methods such as Raman 
and neutron scattering [^. 

For the case of diffusion in the presence of traps distributed randomly with no correla- 
tions, it has been proved rigorously by Donsker and Varadhan that the decay of Po{t) 
with time is slower than exponential, indicating that, in the long time limit, the properties 
of the diffusing particle is dominated by the presence of large trap free regions which have a 
finite (though small) probability of occurence. This is because, if the traps were uniformly 
distributed, the diffusing particle would be trapped at a constant rate, which would result 
in an exponential decay of Po{t). It also means that the quenched disorder average must be 
carried out since an annealed averaging would smear out the trap positions. 

In this paper we ask the question of whether and how this interesting behavior of the 
diffusing particle with traps is modified when we introduce long range correlations in their 
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positions 0. In the equivalent scalar elasticity problem, the correlated traps map to the 
clamping of sites with a correlated distribution. In this sense, the present problem extends 
the so-called fractino problem |^ where a fractal boundary of an otherwise nonfractal object 
is clamped to the case where the bulk of the substrate is itself a fractal. It is also an 
extension of the fracton problem |^ as the traps (or clamped boundaries) are introduced 
into the scalar elasticity of fractals. 

Calculational difficulties have prevented this problem from receiving its share of atten- 
tion even though, often in real physical situations, the positions of traps are correlated at 
least within limited length scales. For the case of uncorrelated trap distribution, a Poisson 
distribution is usually assumed, which simplifies the theoretical calculation P ; however this 
assumption fails for correlated traps. As far as computational calculations, techniques based 
on exact enumeration are highly computer time and memory intensive because of the high 
sensitivity of the behavior of the diffusing particle to the actual positions of the traps in 
the sample, leading to large fluctuations in the measured quantity from sample to sample. 
This translates to requiring an ensemble average to be taken over a large number of disorder 
configurations for meaningful results. Moreover the onset of the asymptotic regime is very 
slow in this type of problem, which adds to the computational difficulties 0. 

In this paper we establish that Po{t) has a stretched exponential form qualitatively similar 
to the case of the Poissonian trap distribution but with a substantially different stretch 
exponent. We show this behavior from the analysis of the density of eigenvalues of the 
transition probability matrix W which describes the diffusion process. We also study the 
finite size scaling of the largest eigenvalue of W, which serves as an extremely powerful tool to 
extract the stretch exponent characterizing the long time behavior of Po{t). A preliminary 
account of this work was presented at the Hayashibara Forum (1995) and appears in its 
proceedings. ||T0| 
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II. IDEAL CHAIN IN CRITICALLY CORRELATED DISORDER 



We choose percolation cluster [11] formed at critical probability of occupation as 
the substrate for diffusion. The sites on the external perimeter (hull) as well as on 
the internal perimeter of the percolation cluster are made absorbing and once the diffusing 
particle reaches the perimeter sites it is trapped permanently. The hull sites and the sites 
constituting the internal perimeter form fractals at pc JT^- Thus placing the traps along the 
perimeters induces spatial correlations among the traps because of the long range correlations 
among the sites of a fractal. While we do not pursue the distinction of external versus internal 
perimeters in this work, there are interesting effects when different boundary conditions are 
applied to them, at least in two dimensions IQ. 

The time evolution of the probability Pi{t) that the diffusing particle is at site i at time t 
is Markovian and the process can be described in the continious time limit by the following 
master equation 

{d/dm{t) = Y.yo,,p,{t) - p,{t), (1) 

i 

where Wij is the hopping rate from site j to i. In this problem we consider only nearest 
neighbour hopping and so Wij is non-zero only for the nearest neighbour pairs i, j. We have 
chosen Wij to be 1/z for all occupied nearest neighbours j, where z is the full coordination 
number of the underlying lattice; otherwise Wij is set equal to zero. What this amounts to is 
that once the particle hops to a trap site there is no further time evolution of that particular 
random walk and only those walks which have escaped getting trapped in the perimeter 
sites evolve further. Thus the traps act as permanent particle absorbers. 

Processes whose time evolution is governed by Eq. (|l|) can be cast into an eigenvalue 
problem of the transition probability matrix W. The matrix elements Wij of W control 
the dynamics of the random walk and the locations of the nonzero elements of W have 
the information about the structure of the underlying fractal substrate responsible for the 
correlations among the trap sites on its perimeters. 
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The density of normal modes of W is related to the return to the starting point prob- 
ability of the diffusing particle by the Laplace transform ||T5|. The probability distribution 
of the number of walks Co{t) which return to their starting point after time t, denoted as 
P{Co{t);t), was studied in and found to be a truncated log-normal distribution. Thus 

P(Co(t);t) ^ ^ exp( —- ). (2) 

In [|l6l it was further found from the first moment of this distribution that 

InPo(t) ~ -t'^'""^"^ (3) 



in the asymptotic long time limit, where Po{t) is obtained as Co{t)/z\ the bar above the 
quantity indicating the quenched disorder average. The exponent Xo is the same one as that 
which appears in the long time behavior of the width of the log-normal distribution, cij^. 
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Since Po{t) has a stretched exponential behavior according to Eq. (|^), we expect the 
density of normal modes p (which is the inverse Laplace transform of Po{t)) to also have a 
stretched exponential form ]TB| , 



Inp(e) ~ -e-'^"/' (5) 
in the limit of small e = | lnA|, where A denotes the eigenvalues of W, and the exponents 



Xo and do are related to each other by do = 4(1 — xo)/i'^Xo ~ !)• The results from [jT6[ gave 
support for the presence of such a behavior although their numerical estimates of do need 
to be improved as they did not take into account the proper normalization of p as well as 
the substantial non-asymptotic effects . 



We can also relate the behavior of p(e) to the finite size effects of the edge of the spectrum. 
First, note that, for any finite substrate however large, the largest eigenvalue of W, which we 
denote by Ai, must be less than one. This is because the eigenvalue of one would indicate the 
existence of a stationary state whereas there are always particles which are trapped, leading 



to a leakage of probability for the diffusing particle in contrast to the diffusion without traps. 
The value of | lnAi|^^ then corresponds to the slowest time scale of the problem, and there 
is always a gap between Ai and one, the latter being what the maximum eigenvalue would 
be for diffusion without traps. The interesting question is whether this gap is bounded from 
below by a nonzero constant or it approaches zero as the substrate size increases. 

Our argument for the behavior Ai — 1 as the substrate size S* — oo follows the observa- 
tion of m for the case of uncorrelated trap distribution. In that case, the return probability 
Po{t) was shown to be a stretched exponential which is slower than a pure exponential (al- 
though not as slow as a power law which would be the case in the absence of traps), and 
this behavior was attributed to the dominance of the trap-free regions although they occur 
relatively infrequently. In our case of critically correlated trap distribution, we also have a 
slow, stretched exponential for Po(t) which we interpret as a similar dominance of the nearly 
trap-free regions. Although there are more traps in the critically correlated case, there are 
also much greater fluctuations in their density. Thus regions of relatively low concentrations 
of the traps are likely to be present, probably with a hierachical size distribution. Since 
the leakage of probability is nearly zero in a nearly trap free region, we would expect Ai to 
approach one as the size of the largest trap free region grows indefinitely. 

Another argument is as follows: if there were indeed a nonzero lower bound for ei = 
I In All, then it would also give a bound for the slowest relaxation rate. Thus Po{t) would 
have to decay at least exponentially in time corresponding to this rate. Since Po{t) in fact 
decays as a stretched exponential with the stretch exponent less than one ( [|1^] as well as 
this work), no such bound can exist for the rate; rather ei — ^ as 5 oo. 

We further propose a scahng relation between ei and S via 



where a possible numerical prefactor has been neglected. This relation follows from the 
assumption that, say, the largest (Ai) and second largest (A2) eigenvalues scale relative to 
the value one in the same way when S ^ 00. That is, if we assume 




(6) 
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r p{e)d{e)^CJ{S) 
Jo 

r p{e)d{e)^C2f{S) 
Jo 



(8) 



(7) 



where 62 = | In A2 1 , then the difference must also scale in the same way, 



/ p(6)rf(6)~(C2-Ci)/(5) 



(9) 



but this latter integral must scale as l/S* if the integral of the density of normal modes p(e) 
is normalized to unity. This means f{S) ~ l/S", thus Eq. (|D follows. 

The assumption of the same asymptotic behavior for the two largest eigenvalues is plau- 
sible if we consider the difference between Ai (or A2) and one to be the reflection of the 
finite size of the largest (or the second largest) trap free region (respectively). As long as 
the large trap free regions are geometrically similar (as would be the case in a hierarchical 
distribution of such regions), and as long as their sizes all go to infinity as S* ^ cxd, it seems 
reasonable that the gaps between Aj {i = 1, 2) and one behave in the same manner. 

If we substitute a stretched exponential form of p(e) in Eq. (^) we get an incomplete 
gamma function. On retaining only the leading term of the incomplete gamma function and 
taking natural log of both sides we get. 



has been assumed and b = c + ln(aa;). 

While the above gives a particular scaling prediction for the case of the critically corre- 
lated trap distribution, the argument leading to Eq. (P) applies more generally. For example, 
a similar procedure should apply to the uncorrelated trap distribution of 0]. Also, for the 



is a stationary state), a similar argument should work with, say, the second (A2) and third 
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(10) 



where a general stretched exponential form of p(e) 




(11) 



trapless case (or the ants [|T5|)) where the largest eigenvalue is actually one (because there 



largest (A3) eigenvalues. Of course, if there is no trap, a probability leakage, per se, does 
not occur. However, each mode with a large A tends to be associated with a blob of high 
connectivity region, to which the same kind of argument can be applied. Indeed, in the 
trap less case, Eq. (|D together with the power law density of states 

p{e) ~ e'^/'-\ (12) 

where dg is the spectral dimension of the substrate, leads to a power law relation between 
S and €2, 

£2-5-2/'^% (13) 
which is identical to the relation proposed on the basis of the finite size scaling of the largest 



nontrivial normal mode and numerically verified in |jT5|. Such finite size scaling relations 
provide a very powerful technique to obtain the quantitative characterization of the density 
of states as they reduce the computational effort drastically, requiring only the information 
on the highest (or second highest) mode. 



III. RESULTS OF NORMAL MODE ANALYSIS 

In this section we give numerical results for the exponent x for the stretched exponential 
decay of the density of states as discussed above. We extract this exponent directly from 
the density of states (cf. Eq. (0)) and also independently from the finite size scaling of the 
largest eigenvalue of W (cf. Eq. (|10|)) in two and three dimensions (square and simple cubic 
lattice, respectively). 

In order to reduce the computational time and memory requirements, we take advantage 
of the fact that we are interested only in the asymptotic long time relaxation of the system 
which is controlled by those normal modes with large eigenvalues A. Thus we use the Arnoldi- 
Saad algorithm |1^,[15| to reduce the original W into a smaller matrix which contains the 
approximate information about the highest normal modes. 
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We analyze the density of states per site, p(e). It is obtained from the eigenvalues of W 
by binning them linearly in the e space. The number of eigenvalues in each bin is divided 
by the bin width, the size of the cluster and also by the number of clusters over which the 
quenched disorder average is performed. The number of independent cluster realizations 
over which the disorder average was taken is of the order 1000 for cluster sizes S = 8000, 
10000, and 50000 and of the order 10000 for S = 5000. We retained in the final results only 
those bins which contained eigenvalues contributed from every substrate realization to avoid 
partial binning. 

Since p(e) has a stretched exponential decay we have plotted — In p(e) versus e in Fig. [l| 
which is expected to have the form ae~^ + c. There is an excellent data collapse for all 
the cluster sizes both in 2d and 3d. The solid curve drawn through the data is obtained 
by fitting the data to an expression of this form. The numerical estimates of x and the 
numerical values of a and c corresponding to the central values of x are tabulated in Table |. 
Note that the value of x cannot be accurately obtained simply from the slope of the double 
log plot of In p(e) with respect to e as the value of the constant c might be appreciable. 

The size of the symbols in Fig. |I] is larger than the cluster to cluster statistical fluctu- 
ations. The quoted error bars are obtained visually from changing the effective estimate 
of X for the nonlinear fit until it no longer fits the data points. The small scaling regime 
which is typical for these kinds of problems (i.e., diffusion with traps) 0] makes the precise 
extraction of x difficult. This accounts for the large error bars in our results as compared to 



the case of diffusion in a percolation cluster without absorbing sites [|T5|. In the latter case 
the density of states has a power law form and the scaling regime increases appreciably with 
the size of the cluster, which resulted in much smaller error bars for the extracted exponents 
even with the same numerical technique. 

In Fig. we plot In 5* with respect to ei in 2d and 3d where 5* is the size of the substrate. 
The solid curves are obtained by fitting the data to an expression aei~^ — (x + 1) Inei + b, 
which is what we expect from the finite size scaling analysis of the largest eigenvalue. The 
size of the symbols are larger than the cluster to cluster fluctuations. The number of clusters 
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over which the disorder average was performed is of the order 1000 for the larger clusters, 
same as for the density of states, and for smaller clusters of size 100, 400, 1000 and 5000 
the average was taken over 10000 clusters. The estimates of x, and the values of a, and b 
corresponding to the central values of x which we obtain from the nonlinear fit are tabulated 
in Table y. The error bar for the value of x is obtained in a similar way to that for the 
density of states. 

We note that the estimates of x and a obtained from the density of states and from 
the finite size scaling of the largest eigenvalue, given in Table | and ||, respectively, are in 
good agreement with each other. The values of do thus obtained, however, turn out to be 
about a factor of two larger than the corresponding stretch exponents for the uncorrelated 
distribution of traps The estimates of x also differ from those of do/2 as given by flG 



due, we believe, primarily to the failure of |jT6[to take proper account of e being not quite in 
the asymptotic region (thus proper normalization and prefactors becoming important) (cf. 
p!7|). However, the exponent Xo (Eq. (^) is relatively insensitive to x (= do/2) and the 
analysis of Po{t) in |]TB|] is not affected by these problems. We also believe that their main 
conclusions remain valid. 



IV. SUMMARY 

In summary, we have studied the problem of diffusion on a substrate with permanent 
traps which are distributed with critical correlation in their positions. The critical correlation 
has been modeled by placing the traps on the perimeters of critical percolation clusters, 
which are obtained by Monte Carlo simulation. The statistical behavior of the diffusing 
particle in the time domain is then mapped to the scalar elasticity problem with fixed, 
fractal boundaries, and a numerical analysis of the vibrational density of states of the latter 
problem is carried out using an approximate diagonalization algorithm of Arnoldi and Saad 
p!8| . This problem may be considered a generalization both of the fracton problem of 
Alexander and Orbach where there are no traps but the substrate is fractal, and of the 
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fractino problem of Sapoval et al where a fractal boundary is clamped (equivalent to 
traps) but the substrate itself is not fractal. 

We have shown that introducing critical spatial correlations among the traps does not 
change the qualitative behavior of the density of states as compared to the case of uncorre- 
lated traps, but that there are substantial effects of correlations quantitatively. This strongly 
suggests that the diffusion process is dictated by the presence of nearly trapless regions sim- 
ilarly to the case of uncorrelated traps, but the long range correlations among the trap 
positions lead to a decrease in the number and size of these regions nearly free of traps, 
which induces a faster decay in the density of states and consequently a smaller probability 
of return to the starting point after time t. This is reflected in a much larger exponent in 
the stretched exponential form of density of states than the case of uncorrelated traps. 

On the other hand, as compared to the case of the same fractal substrate but with no 
traps, the density of states is qualitatively different, since the latter problem produces a 
power law density of states in the low energy limit known as fractons, accumulating to an 
infinite density toward the maximum eigenvalue of one ||TB[. In contrast, with the perimeters 
acting as traps, the density of states becomes a stretched exponential, rapidly falling to zero 
toward the maximum eigenvalue. 
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FIGURES 

FIG. 1. Natural log of the density of eigenvalues p{e) against e in d = 2 and d = 3. The data 
collapse very well for different sizes of substrates. 

FIG. 2. Natural log of the size of substrate \nS against ei in d = 2 and d = 3. 
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TABLES 

TABLE I. Estimates of the exponent x and the constants a and c from p{e) for the square 
lattice in two dimensions and simple cubic lattice in three dimensions. The data for — In p(e) have 
been fitted to the form ae~^ + c. 
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TABLE IL Estimates of the exponent x and the constants a and h = In(ax) ± c from the 
finite size scaling of ei, for the square lattice in two dimensions and simple cubic lattice in three 
dimensions. The data have been fitted to the form ae"^ — (x + 1) In e + 6. 
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